library(sasld)

# -------------------------------------------------------------------
# ATTENZIONE: la funzione sample e' stata modificata 
# a partire da R 3.6.0 per cui i "vecchi risultati"
# non posso essere riprodotti di default
# http://hostingwin.unitn.it/micciolo/PeM/01-rng.pdf
# -------------------------------------------------------------------
# se si vogliono riprodurre i risultati dei Laboratori
# con le versioni pił recenti di R va eseguita la seguente istruzione
# -------------------------------------------------------------------
RNGkind(sample.kind="Rounding")
# -------------------------------------------------------------------


# 6.1 ---------------------------------------------------------------
tl <- c("R","S")

set.seed(365274)
paste(sample(tl,20,replace=TRUE,prob=c(0.5,0.5)),collapse=" ")

# facciamolo 10 volte
set.seed(365274)
for (i in 1:10) {
    out <- sample(tl,20,replace=TRUE,prob=c(0.5,0.5))

    x <- (out == "R")*1
    cat(paste(out,collapse=" ")," - ",shots(x),"\n")
}



ris <- numeric(10000)
set.seed(365274)
for (i in 1:10000) {
    out <- sample(tl,20,replace=TRUE,prob=c(0.5,0.5))
    x <- (out == "R")*1
    ris[i] <- shots(x)
}
table(ris)
round(table(ris)/10000,2)

tbl <- table(ris)
sum(tbl[10:15])

0.5^20

mean(ris)


a <- c("A","B")
set.seed(123456)
x <- 0
nA <- 0
nB <- 0
ris <- character(7)
repeat {
    y <- sample(a,1)
    x <- x+1
    ris[x] <- y
    ifelse((y == "A"),nA <- nA+1,nB <- nB + 1)
    if((nA == 4) | (nB == 4)) break
}
ris
x


out <- numeric(10000)
set.seed(123456)
for (i in 1:10000) {
    x <- 0
    nA <- 0
    nB <- 0
    ris <- character(7)
    repeat
{
        y <- sample(a,1)
        x <- x+1
        ris[x] <- y
        ifelse((y == "A"),nA <- nA+1,nB <- nB + 1)
        if((nA == 4) | (nB == 4)) break
}
    out[i] <- x
}

table(out)

tbl <- table(out)
(4*tbl[1] + 5*tbl[2] + 6*tbl[3] + 7*tbl[4])/10000

0.5^4*2


out <- numeric(10000)
set.seed(654321)
for (i in 1:10000) {
x <- 0
nA <- 0
nB <- 0
ris <- character(7)
repeat
{
    y <- sample(a,1,prob=c(0.9,0.1))
    x <- x+1
    ris[x] <- y
    ifelse((y == "A"),nA <- nA+1,nB <- nB + 1)
    if((nA == 4) | (nB == 4)) break
}
out[i] <- x
}
table(out)

tbl <- table(out)
(4*tbl[1] + 5*tbl[2] + 6*tbl[3] + 7*tbl[4])/10000

0.9^4 + 0.1^4
# -------------------------------------------------------------------




# 6.2 Normale -------------------------------------------------------
curve(dnorm(x,mean=65,sd=3.5),55,85)
curve(dnorm(x,mean=70,sd=4.0),55,85,add=TRUE)

Normale(mu=65, sigma=3.5)
Normale(mu=70, sigma=4.0, add=TRUE)

x <- seq(1.4,1.49,0.01)
pnorm(x)

pnorm(60,mean=65,sd=3.5)
pnorm(70,mean=65,sd=3.5)
pnorm(70,mean=65,sd=3.5) - pnorm(60,mean=65,sd=3.5)

ProbNorm(da=60,a=70,mu=65,sigma=3.5,col="grey")
(60-65)/3.5
(70-65)/3.5

ProbNorm(da=80,a=90,mu=83,sigma=5)

1-pnorm(650,mean=500,sd=100)
pnorm(650,mean=500,sd=100,lower.tail=FALSE)

qnorm(0.98)
qnorm(0.98,mean=100,sd=16)
# -------------------------------------------------------------------




# 6.3 ---------------------------------------------------------------
dbinom(c(0,1,2,3),size=3,prob=0.2)

dbinom(c(0:10),size=10,prob=0.5)

Binomiale(n=10,p=0.5)

Binomiale(n=10,p=0.5,da=0,a=0,grafico=FALSE)

sum(dbinom(c(444500:500000),size=500000,p=0.554))
sum(dbinom(c(444500:500000),size=500000,p=0.554,log=TRUE))
sum(dbinom(c(444500:500000),size=500000,p=0.554,log=TRUE))*log10(exp(1))

m <- 500000*0.554
v <- 500000*0.554*0.446
z <- (444500-m)/sqrt(v)
pnorm(476.5,lower.tail=FALSE)
pnorm(476.5,lower.tail=FALSE,log=TRUE)
pnorm(476.5,lower.tail=FALSE,log=TRUE)*log10(exp(1))
# -------------------------------------------------------------------




# 6.4 ---------------------------------------------------------------
dbinom(81,162,0.42)

sum(dbinom(c(81:162),162,0.42))

mu <- 162*0.42
va <- 162*0.42*0.58

pnorm(81,68.04,sqrt(39.4632),lower.tail=FALSE)

pnorm(80.5,68.04,sqrt(39.4632),lower.tail=FALSE)

p81 <- sum(dbinom(c(81:162),162,0.42))
x <- c(0:1000)
p <- p81*(1-p81)^x

x <- c(0:1000)
p <- 0.024*0.976^x

sum(p)

sum(x*p)

(1-0.024)/0.024
# -------------------------------------------------------------------